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We consider the stability and dynamics of multiple dark solitons in cigar-shaped Bose-Einstein condensates 
(BECs). Our study is motivated by the fact that multiple matter-wave dark solitons may naturally form in such 
settings as per our recent work [Phys. Rev. Lett. 101, 130401 (2008)]. First, we study the dark soliton inter- 
actions and show that the dynamics of well-separated solitons (i.e., ones that undergo a collision with relatively 
low velocities) can be analyzed by means of particle-like equations of motion. The latter take into regard the 
repulsion between solitons (via an effective repulsive potential) and the confinement and dimensionality of the 
system (via an effective parabolic trap for each soliton). Next, based on the fact that stationary, well-separated 
dark multi-soliton states emerge as a nonlinear continuation of the appropriate excited eigensates of the quan- 
tum harmonic oscillator, we use a Bogoliubov-de Gennes analysis to systematically study the stability of such 
structures. We find that for a sufficiently large number of atoms, multiple soliton states may be dynamically 
stable, while for a small number of atoms, we predict a dynamical instability emerging from resonance effects 
between the eigenfrequencies of the soliton modes and the intrinsic excitation frequencies of the condensate. 
Finally we present experimental realizations of multi-soliton states including a three-soliton state consisting of 
two solitons oscillating around a stationary one. 

PACS numbers: 



I. INTRODUCTION 



Dark solitons, namely localized density dips on top of a stable continuous-wave background (or a background of finite extent) 
with a phase jump across their density minimum, are fundamental envelope solitons supported in nonlinear dispersive media. 
These nonlinear waves emerge in media with a positive (negative) group-velocity dispersion and defocusing (focusing) nonlin- 
earity, with a proper model describing their evolution being the so-called defocusing nonlinear Schrodinger (NLS) equation pi- 
Dark solitons have been studied extensively in the field of nonlinear optics (from which the term "dark" was coined) S, but 
they have also been observed in diverse physical contexts, including liquids H, mechanical systems [4], magnetic films j5|], and 
so on. 

More recently, dark solitons have attracted much attention in the physics of Bose-Einstein condensates (BECs) 06 ,01, where 
they appear as fundamental macroscopic nonlinear excitations of BECs with repulsive interatomic interactions K , M l. It is, 
therefore, not surprising that experimental work on matter-wave dark solitons started as early as ten years ago ITolfTll L31illH 
and, very recently, continued even more intensively 11151 llo, Il7l Il8l Il9ll . These important "new age" experiments highlighted 
various salient features of dark solitons, verified previous theoretical predictions and offered motivation for further investigations. 
A pertinent example is the creation of more than one matter-wave dark solitons 03 H (see also Ref. 01) in cigar-shaped 
condensates, which were allowed to interact. This invites a revisiting of the topic of dark solitons and especially of their 
interactions in the particular context of BECs; the latter has a number of particularities including e.g., the confinement that is 
routinely used to trap and cool the atomic cloud (a, |7[] . 

Multiple dark soliton solutions of the defocusing NLS equation were first obtained in Ref. 01 by means of the inverse 
scattering method. Later, an analytical form of a solution of the NLS equation composed of two dark solitons of different depths 
and velocities was found ll20ll (see also the more recent works 1I2H l22ll ). and it was shown that the interaction between dark 
solitons is repulsive. Subsequent theoretical studies focused on the interactions and collisions of dark solitons in the context 
of nonlinear optics I23I l24l l25ll and later in BECs lf26ll . while relevant experimental results (see Refs. ll27ll for optical dark 
solitons and P 4, 1 1711 for atomic dark solitons) also examined the interaction between two dark solitons. However, following the 
recent experimental methodology of Ref. lfl8ll . it is in principle possible to generate multiple (in fact, in principle, an arbitrary 
number of) dark solitons: this can be done by releasing a BEC from a double-well potential into a harmonic trap within the 
experimentally accessible, so-called, dimensionality crossover regime between one-dimension (ID) and three-dimensions (3D) 
Pal . In such a case, it is clear that a study of multiple matter- wave dark solitons, and their interactions, should be performed in 
a theoretical framework that takes into regard basic features of the pertinent experiment, such as the effect of dimensionality and 
the corresponding effective confinement of the condensate. 

In this work, our scope is to analyze this problem, namely the statics and dynamics of multiple matter-wave dark solitons 
in cigar-shaped condensates. Based on the fact that recent atomic dark soliton experiments were performed at extremely low 
temperatures and with sufficiently large number of atoms, we may safely adopt a mean-field theoretical approach. In particular, 
we will perform our analysis in the framework of the effectively ID Gross-Pitaevskii (GP) equation with a non-cubic nonlinearity 
that was first presented in Ref. [28] and later was also derived and tested in Refs. B29I1 [this is a distinguishing feature of our work 
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with respect to most of the above references which considered dark soliton interactions in the standard homogeneous defocusing 
cubic NLS setting]. Our analysis starts by considering the dynamics of multiple dark solitons which is studied as follows. First, 
we consider the weakly-interacting limit of the non-cubic GP equation (namely the traditional defocusing NLS model) and, in 
the absence of the trap, we derive an effective repulsive potential for the interaction between two solitons. It is shown that this 
potential can successfully be used to describe the interactions between "low-speed" solitons (with velocities less than the half of 
the speed of sound). Such solitons are, "well-separated" in the sense that they always can be identified as distinguishable objects, 
even at the collision point. Then, using this potential, we obtain a set of equations of motion for the coordinates of an arbitrary 
number of solitons. Our approach, is finally applied to the full problem under consideration (with the non-cubic nonlinearity and 
the external harmonic trap), upon incorporating an effective harmonic potential with a corresponding characteristic frequency: 
this is actually the eigenfrequency of the first anomalous mode of the system 0], corresponding to the oscillation frequency 
of a single dark soliton in the trap (see relevant results in Refs. QOll and Bill for the purely-lD and dimensionality-crossover 
regimes). Such an ad hoc decomposition of the principal physical mechanisms affecting the solitons was first introduced in Ref. 
01 (for the case of two symmetrically interacting dark solitons), and will be validated a posteriori herein, by means of direct 
numerical simulations. 

The above methodology for the study of multiple atomic dark solitons is directly connected to the Bogoliubov-de Gennes 
(BdG) spectrum of excitations of stationary dark soliton states. The latter are obtained when linearizing around the nonlinear 
counterparts of the respective linear states (corresponding to the eigenmodes of the quantum harmonic oscillator) ll32ll and their 
properties are studied by means of the well-known BdG equations 0]. Such an analysis reveals that the spectrum of the rt-th 
excited state consists of one zero eigenvalue, n double eigenvalues (accounted for by the presence of the harmonic trap), and 
infinitely many simple ones. In the nonlinear regime, one of the eigenvalues of each double pair possesses a topological property 
of, so-called, negative energy (in the physical literature) lf33tl or negative Krein signature (in the mathematical literature) HH; 
practically, this means that it becomes structurally unstable, i.e., it becomes complex, upon collision with other eigenvalues. The 
eigenvalues with negative Krein signature are actually associated with the anomalous modes [7] appearing in the BdG spectrum. 
In our case of multiple dark solitons, the number of anomalous modes in the excitation spectrum equals to the number of dark 
solitons ll35ll . which is in agreement with the fact that the number of eigenvalues with negative Krein signature equals to the 
number of the nodes of the stationary state [36]. More generally, we conjecture (based on the results below for the cases of 
two- and three-solitons) that in the case of an rc-dark soliton sequence (pertinent to an n-th order nonlinear state), the anomalous 
modes of the system correspond to the excitation of the normal modes of the "dark-soliton lattice". 

The paper is organized as follows. In section II we present the model and make some general remarks on the theoretical setup. 
Section III is devoted to the dynamics of multiple solitons. In particular, first we analyze the homogeneous weakly-interacting 
case, and derive the effective repulsive potential for two solitons undergoing a symmetric collision. Then, we generalize these 
results to include the cases of asymmetric collisions and multiple solitons, as well as to tackle the full problem, taking into regard 
the external harmonic trap and the dimensionality of the condensate. In section IV we study the stability of the stationary multi- 
soliton states via a BdG analysis. We analyze, in particular, the pertinent Bogoliubov spectra, paying special attention to the 
anomalous modes of the system. We illustrate how these anomalous modes correspond to "normal modes" of the "dark-soliton- 
lattice", e.g., in-phase and out-of-phase oscillating dark soliton states. We also predict the onset of dynamical instabilities due to 
resonance between the eigenfrequencies of these normal modes and the excitation frequencies of the background condensate. We 
illustrate under what conditions such instabilities may be observed in future experiments. In Section V we present experimental 
realizations of multi-soliton states including a three-soliton state consisting of two solitons oscillating around a stationary one. 
Section VI concludes the paper, summarizing our findings and presenting some directions of future study. 



II. THE MODEL AND THEORETICAL SETUP 



We consider a BEC confined in a highly elongated trap, with longitudinal and transverse confining frequencies (denoted by 
uj z and lo±, respectively) such that uj z <C uj±. In this case, it can be found ll28ll29ll that use of the adiabatic approximation, in 
combination with a variational approach for determining the local transverse chemical potential, leads to the following effectively 
ID GP equation, 



di[> 



Ira dz 2 



+ V ext {z) + hw^y/l + 4a|i/f 



1>, 



(1) 



where i/j(z, t) is the longitudinal part of the condensate's wave function normalized to the number of atoms, i.e., N = 
f_°° \ip\ 2 dx, a is the s-wave scattering length, m is the atomic mass, and V ex t{z) is the longitudinal part of the external 
trapping potential, assumed to take the standard harmonic form V ex t(z) — (l/2)mui 2 z 2 . As demonstrated in Refs. H29fl . Eq. 
(fl]l provides accurate results in the dimensionality crossover and the Thomas-Fermi limit, thus describing the axial dynamics of 
cigar-shaped BECs in a very good approximation to the 3D Gross-Pitaevskii (GP) equation. Notice that in the weakly-interacting 
limit, 4a|?/'| 2 <C 1, Eq. (Q]) is reduced to the usual ID GP equation with a cubic nonlinearity, characterized by an effective ID 
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coupling constant guj = 2ahui±. Equation (0 can be expressed in the following dimensionless form, 



.dip 



1 d 2 1 



2dz2 ■ 2 -^ 2 + VTTw 



^, (2) 

where = lu z /lu± is the normalized trap strength. In Eq. (0, the density \ijj\ 2 is measured in units of 2aN, while length, time 
and energy are measured, respectively, in units of the transverse harmonic oscillator length a± = y / h/mu!±, w7 , and fruj±. 

Exact analytical dark soliton solutions of Eq. are not available. However, following the lines of Ref. 113711 where a NLS 
equation with a generalized defocusing nonlinearity was considered, dark soliton solutions can be found in an implicit form (via 
a phase-plane analysis) or in an approximate form (via the small-amplitude approximation). On the other hand, exact analytical 
dark soliton solutions are available in the above mentioned weakly-interacting limit (2|?/>| 2 <C 1), and in the absence of the 
external potential, i.e., for the cubically nonlinear defocusing NLS model. A single dark soliton solution on top of a background 
with constant density n = no = fJb (with ji being the chemical potential) has the form QJJ] , 

ip(z, t) = y/no \iv + B tanh(?7)] exp(— i/it), (3) 

where r\ = ^fn^B(z — y/novt), the parameter B = \/l — v 2 sets the soliton depth, ^/tlqB, while the parameter v sets the 
soliton velocity, y/rujv. Note that for v = the dark soliton becomes a stationary kink (alias "black" soliton), while for v = 1 
the dark soliton solution ([3]) becomes the background solution. Multiple dark soliton solutions are also available GHUH. In 
the simplest case of a two-soliton solution, with the two solitons moving with equal velocities, v\ = — v 2 — v, the wave function 
can be expressed as 112 111 (see also B22I1 ): 

Hz,t) = ^p^-exp(-int), (4) 
G(z, t) 

whereF = 2(tiq — 2n m i n ) cosh(T) — 2nof cosh(Z)+i sinh(T), G = 2 v / no cosh(T) + 2^/n m i„ cosh(Z), while Z = 2^/n^Bz, 
T = 2y / 'n m i n (no — n m i n )t, and n m i n — tiq — jiqB 2 = riQi/ 2 is the minimum density (i.e., the density at the center of each 
soliton). This equation will be useful for the analysis of dark soliton interactions (see next section). 

Generally, the single dark soliton, as well as all higher-order dark soliton states, can be obtained in a stationary form from the 
non-interacting (linear) limit of Eq. (0, corresponding to N — > 0. In this case, Eq. |2]is reduced to a linear Schrodinger equation 
for a confined single-particle state, and this problem becomes the equation for the quantum harmonic oscillator; the latter, is 
characterized by discrete energy levels and corresponding localized eigenmodes described by the Hermite-Gauss polynomials 
dUt]. Then, in the weakly-interacting case, where Eq. (fJJ becomes the cubic NLS equation, all these eigenmodes exist for the 
nonlinear problem as well fl32ll . describing an analytical continuation of the linear modes to a set of nonlinear stationary states. 
Additionally, recent analysis and numerical results [ 39] suggest that there are no solutions of Eq. (0 without a linear counterpart. 
In fact, the effect of interactions (i.e., the effect of nonlinearity) transforms all higher-order stationary modes into a sequence of 
parabolically confined dark solitons I32I1 . From the physical point of view, the higher-order stationary modes exist due to the fact 
that the repulsion between dark solitons |2(A |2U l23l u4. 12511 is counter-balanced by the restoring force induced by the trapping 
potential. 

Below, we are going to analyze the stability of nonlinear modes (namely stationary multi-dark-soliton states) of Eq. (0 by 
means of the BdG equations. In particular, first we identify a numerically exact (up to a prescribed tolerance), stationary soliton 
state, ipT>s(x), using a fixed point algorithm (e.g., a Newton-Raphson method). Then, considering small perturbations of this 
state of the form, 



1>{z, t) = [ifoa(z) + (u(z)e- iut + v*(z)e 1u ' t ) 



e~ vt , (5) 



(where * denotes complex conjugate), we derive from Eq. 0l the following BdG equations: 

[H - fx + f]u + gv = uju, (6) 
[H — fx + f]v + gu = —uv, (7) 

where H = —(l/2)d 2 + (l/2)f2 2 z 2 is the single particle operator, fx is the chemical potential and the functions / and g are 
given by / = g + a/1 + 4tids. an d 9 — 2"0ds/V1 + 4tids (with n^s = IV'ds) 2 )- Then, solving Eqs. ©-(0, we are going to 
find the eigenfrequencies lu = uo r + iu>i and the amplitudes u and v of the normal modes of the system. Note that due to the 
Hamiltonian nature of the system, if oj is an eigenfrequency of the Bogoliubov spectrum, so are —lu, lu* and — lu*. Notice that a 
linearly stable configuration is tantamount to lui = 0, i.e., all eigenfrequencies being real. 

The stability of nonlinear modes has already been considered in several works |40 l l4ll l42ll in the framework of the ID GP 
equation. According to Ref. 1I4()I1 stable nonlinear modes exist, but results of Ref. [42], obtained near the non-interacting limit, 
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appear to contradict those of Ref . l40ll . In fact, in Ref . J42ll it was claimed that, apart from the first one, all higher-order nonlinear 
modes are unstable for repulsive BECs (but can be stabilized by using anharmonic traps). The mechanism of this instability can 
be interpreted as follows: in the non-interacting (linear) limit, the excitation spectrum of the n-th mode possesses n equidistant 
(due to the harmonic trap) double eigenvalues. Then, passing into the nonlinear interacting regime, the negative energy of one 
member of each of these double eigenfrequency pairs practically leads to structural instability, i.e., each pair and its opposite 
split into a complex quartet of eigenfrequencies. 

Below we will show that the above mentioned instability of higher-order nonlinear modes occurs near the non-interacting 
limit, but will cease to exist sufficiently deep inside the nonlinear regime (i.e., for sufficiently large condensates with large 
numbers of atoms TV). In fact, considering both the second- and the third-order nonlinear mode of Eq. (0, we will demonstrate 
that they are initially (i.e., for a small number of atoms, sufficiently close to the linear limit) linearly unstable. Nevertheless, 
numerical continuation of the corresponding waveforms to larger values of N, reveals a critical value of the number of atoms 
(depending on the anisotropy of the external harmonic trap) which, if exceeded, all the eigenfrequencies become real and, thus, 
the nonlinear stationary state becomes linearly stable. 



III. DYNAMICS OF MULTIPLE DARK SOLITONS 

A. Multiple soliton interactions - The homogeneous case 

Let us consider, at first, the simplest possible multi-soliton state, namely a pair of dark solitons located at z = ±zq, and 
moving with opposite velocities, i.e., v\ = —v-x = v. Then, in the framework of the weakly-interacting limit of Eq. (0 (and in 
the absence of an external potential), we may derive an equation for the trajectory of the soliton coordinate, z , as a function of 
time. This can be done upon identifying the soliton coordinate as the location of the minimum density, and using the equation 
d\ifj\ 2 /dz = 0, with ip(z, t) given in Eq. (01), to obtain the result: 



cosh(2V^Bz ) = J^cosh(T) - 2./ """" } (8) 
• - mn V ri cosh(T) 



(recall that T = 2-^/n m i„(no — n m i n )t). Then, Eq. ^ can be used for the determination of the distance 2zq between the two 
solitons at the point of their closest proximity (corresponding to t = 0): 



' icosbr 1 ( - 2./^^- I . (9) 



This equation shows that the closest proximity distance becomes zero for n min /n n = v 2 = 1/4. Physically, this means that 
there exists a critical value of the soliton velocity, namely v c = 1/2, which separates two different regimes: in the first regime, 
"high-speed" solitons with v > v c are transmitted through each other [their high kinetic energy overcomes the interparticle 
repulsion], while in the second regime, "low-speed" solitons with v < v c are reflected by each other [their lower kinetic energy 
in this regime is insufficient to overcome the interparticle repulsion]. In fact, as seen in the panels of Fig. Q] the wave function 
of the low-speed (high-speed) solitons exhibits two separate minima (a single non-zero minimum) at the collision point, namely 
\ij){zQ,t = 0)| 2 = (\ip(zQ,t = 0)| 2 7^ 0). Note that in the case of the critical velocity v c — 1/2, the two-soliton wave 
function exhibits a single zero minimum at the collision point. The above analysis underscores the fact that low-speed solitons 
are actually "well-separated" solitons, in the sense that they can always be characterized by two individual density minima even 
at the collision point (the point of their closest proximity). On the contrary, the high-speed solitons completely overlap at the 
collision point and, thus, are not distinguishable during the collision. Thus, well-separated solitons appear to be reflected by each 
other and can safely be regarded as hard-sphere-like particles that interact through an effective repulsive potential (although, as 
we will see quantitatively, the description below will be surprisingly accurate even in the non-well-separated case). 

To further elaborate on the above, let us consider the limiting case of extremely slow solitons, namely no /n m in — 1/2 <^ 1 /4, 
for which the soliton separation is large for every time, i.e., the closest proximity distance is Zq 3> 0. In this case, the second 
term in the right-hand side of Eq.® is much smaller than the first one for every time (including t = 0) and can be ignored. This 
way, the soliton coordinate can be expressed as: 

zq = cosh -1 [v 1 cosh(2no^-Bi)l , (10) 

which yields the soliton velocities: 

dz — sinh(2n iABi) 

-TT = V^o I (11) 
<v- x cos\i 2 (2n vBt) - 1 
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FIG. 1: Top panels: Soliton trajectories for symmetric two-soliton collisions for different initial velocities: v — 0.2 (left) and v = 0.8 (right). 
In the former (latter) case the solitons are reflected by (transmitted through) each other. Bottom panels: Soliton trajectories for an asymmetric 
two-soliton collision for initial velocities v\ — 0.5 and 1/2 = (left), and for a three-soliton collision for initial velocities v\ = —^3 = 0.4 
and ^2 = (right). In all panels, shown are density plots of the wave function carrying the solitons as obtained by direct numerical integration 
of the homogeneous NLS equation. The solid lines correspond to the solution of Eq. d 1 4b (i.e., employing the effective interaction potential). 
In all cases the normalized chemical potential is \i — 1. 



The above equation shows that in the limit t — ► ±00, the soliton velocities take the asymptotic values dzg/dt = ±-JnoU, namely 
the values of the velocities of each individual soliton [see the definition of the single soliton velocity beneath Eq. (0]. On the 
other hand, at t = 0, Eq. (fTTb yields dzo/dt = 0; this means that as the dark solitons are approaching each other, they become 
slower, i.e., darker, and at t — (corresponding to the point of their closest proximity) they become black, remaining at some 
distance away from each other. After such a, so-called, head-on "black collision" [26], the dark solitons are reflected by each 
other and continue their motion in opposite directions. 

We now proceed to determine the effective repulsive potential for well-separated dark solitons. This can be done by deter- 
mining, at first, an equation of motion for the soliton coordinate: differentiating Eq. ( TTOb twice with respect to time, and using 
Eq. ([8]) (without the second term, which is negligible for well-separated solitons), we obtain the equation of motion in the form 
d 2 zo/dt 2 = —dV(zo)/dzo, with the repulsive potential being given by: 

V(z ) = l ^ . (12) 

2 sinh (2y/rioBzo) 

It is worth noting here that since B = y/1 — v 2 , the above potential is, in principle, a velocity dependent one. Note that Eq. (fl~2b 
recovers the result obtained in Ref. l25ll by means of a Lagrangian approach (in that work, the sinh term in the denominator 
appears as a cosh term due to a typographical error ll43to . 

Although the potential of Eq. ( fT2l is formally applicable only to symmetric collisions, it can nevertheless be applied also in 
the case of non-symmetric collisions provided that an "average depth" of the two solitons is employed. In fact, it is possible 
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to generalize this concept for an arbitrary number of soli tons, n: assuming that the i-th soliton (with i = 1,2, ••• , n) is 
characterized by a darkness Bi, velocity = wl — B 2 , and a position z. ; , we may define the average depth B t j = (l/2)(S i + 
Bj) and the relative coordinate = (1/2) (zi — Zj) for solitons i and j, and express the interaction potential Vi in the presence 
of other solitons, as: 

v l = y — 2 noB - . as) 

2smh [y/n^Bijizi - Zj)\ 

Notice that Eq. ( fT3l is reduced to Eq. ( fT2T > for Vi = — fj = v (i.e., for Bi = Bj = B) and Zi — Zj = 2zo. 

Using Eq. ( fT3l l, it is now straightforward to obtain equations of motion for a "lattice" consisting of an arbitrary number of 
dark solitons. Taking into regard that the Lagrangian L of n interacting solitons is L = T — V, where T = Y17=i (l/2)ii 2 (with 
ii = dzi/dt) and V = J27=i 316 tne kinetic and potential energy, respectively, the Euler-Lagrange equations, d(di i L)/dt — 
d Zi L = 0, lead to the following set of dynamical evolution equations: 

d 2 V . d 2 V ,.\ dV n 

Zk + 7T. n- Zk I ~n = (14) 



f— ' \ozkdzi ozkdzi J dz., 

k—l 

These n coupled equations of motion can then be used to calculate the trajectories zt(t) of n-interacting dark solitons. It is worth 
pointing out here that in deriving Eqs. ( TBI , we have attempted to incorporate the character of the solitary waves as "deformable 
particles" with a velocity-dependent interaction potential (cf. Eq. ( fT2b and related discussion). This approximation will be tested 
a posteriori through the detailed comparison of the particle-based and the GP-based dynamical results. 

We have performed systematic numerical simulations to investigate the range of validity of Eqs. (TBI) , both for cases of 
symmetric and asymmetric soliton colisions, as well as both for cases of low-speed (well-separated) and high-speed dark solitons. 
Various relevant examples are shown in Fig. Q] The simulations confirm that as long as the dark solitons are well-separated from 
each other, i.e., if their depth (velocities) is (are) sufficiently large (small), their trajectories found by means of Eqs. ( TPfl ) almost 
coincide with the ones found by direct numerical integration of the NLS equation. This excellent agreement can be illustrated not 
only qualitatively but also quantitatively: this can be done upon comparing the exact results for the collision-induced phase-shifts 
of the soliton trajectories to the ones found numerically by means of Eq. ( TBI . In the case of two solitons, these phase-shifts 
were calculated analytically in Ref . [ 1 ] and have the following form, 

= ^ { r- V2 l +[Bi+B2 L 

2B 1 (isi-V2) 2 + (B 1 -B 2 y 

s Z2 = -jL ta (*-»»);+(*+*);. ( i6) 

2B 2 {vi - v 2 ? + (Si - B 2 ) 2 

In Fig. |2]we compare the exact phase-shifts provided by the above expressions to the ones determined by means of Eq. ( TBI ), 
which employ the effective repulsive potential of Eq. (fT2] i. Both cases of a symmetric (top panels) and an asymmetric (bottom 
left panel) collision are shown; it is clearly observed that the agreement between the two approaches is very good for soliton 
velocities v < 1/2, i.e., for well-separated dark solitons. Notice that in the case of an asymmetric collision, v\ ^ v 2 , the two 
shifts are not equal, \&z\\ ^ \5z 2 \, while in the case of a symmetric collision, v\ = — v 2 = v (and, thus, B\ = B 2 = B), the 
phase shifts become equal, \8zi\ — \8z 2 \ — (2B)^ 1 ln(l + B 2 /u 2 ). 



B. Dynamics and interactions of multiple solitons in the trap 

The above analysis of soliton interactions in the homogeneous case is of use in the inhomogeneous case as well. In particular, 
here we will consider multiple dark solitons in the presence of an harmonic trap, also taking into regard that the condensate 
is cigar-shaped, so that the proper model is Eq. (f2]l (rather than its weakly-interacting limiting case, i.e., the usual cubic NLS 
equation considered above). In such an experimentally relevant situation, we may employ the theoretical approach adopted in 
Ref. S and use an interaction potential for dark solitons of the form: 

V* ff = vS f (z i )+V i (z i) z i ), (17) 

where Vi(zi, Zi) is the interaction potential of Eq. ( fT3l and the effective trapping potential V^J/ is given by: 

v:i t } {zi) = \u:i sc z 2 . (is) 
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FIG. 2: Collision-induced phase-shift of the solitons as a function of the velocity v. Solid lines depict the analytical result of Eqs. d 1 5b-d 1 6fc 
and the dashed lines show the results obtained from Eq. d!4t . which employs the effective repulsive potential. The agreement between the two 
approaches is very good for well-separated solitons, i.e., for v < 1/2. The left panel shows the case of a symmetric collision, and the right one 
the case of an asymmetric collision of a moving soliton with v = V\ and a stationary soliton (i.e., a black soliton with ^2=0- see bottom 
left panel of Fig. [TJ. The dashed vertical line depicts the critical value v — 1/2 and defines the range of applicability of the approach based on 
the effective repulsive potential. 



The underlying assumption within this decomposition is that dark soliton is an effective particle moving under the combined 
influence of external forces from the confining potential and from the other solitons within the configuration. Each of these 
individual forces has an associated potential and hence the total force and associated motion stem from the combination of these 
potentials. A more subtle assumption is that while the effect of dimensionality on a single soliton is captured in the effective 
oscillation frequency lo osc discussed in detail in the next paragraph, the tail-tail interaction of the waves is well approximated 
by its NLS counterpart. These assumptions will be validated a posteriori through our comparisons between theoretical and 
numerical results below. 

Here, lo osc is the oscillation frequency of a single dark soliton in the harmonic trap, which coincides with the lowest anomalous 
mode of the system [31]. In the purely ID regime, and for sufficiently large numberof atoms (i.e., in the so-called Thomas-Fermi 
regime the soliton oscillation frequency is iv osc = Sl/y/2, where f2 is the normalized trap strength. This can be derived 

either by analyzing the dynamics of the dark soliton in the trap (see, e.g., |@] and references therein) or by means of a BdG 
analysis 13011 . In the case of cigar-shaped BECs under consideration, the oscillation frequency is upshifted [31] and, generally, 
takes values in the interval Q / v2 < lu osc < Q. 

Based on the above discussion, the interaction potential of Eq. ([T8l takes into account both the effective harmonic trap 
(including the dimensionality of the system), V^J/(zi), and the inter-soliton interaction potential, Vi(zi,Zi) (derived for the 
homogeneous ID regime). This potential has already been successfully used in Ref. 01811 . where the experimental findings for 
the symmetric collisions between two dark solitons were found to be in excellent agreement with the corresponding theoretical 
results. Here, we will show that the approach based on the use of the effective potential of Eq. ( TPTI i can also be generalized to 
the case of asymmetric collisions. Such a case can also be investigated experimentally in the context of the experimental setup 
of lfl8ll . In particular, in [18] an even number of dark solitons was created by merging two condensates initially prepared in a 
double-well trap, with a zero-phase difference between them. In principle, it is also easy to create an odd number of solitons, 
upon introducing a non-zero phase difference between the wells, which would lead to an asymmetric evolution pattern of the 
solitary waves. If, furthermore, the phase difference is exactly equal to n, then a stationary (black) dark soliton is created exactly 
at the center of the harmonic trap. In fact, this procedure has already been used in relevant interference experiments (not directly 
connected, however, to dark solitons) ll45ll . 

In Fig. [3] we show an example of an asymmetric collision between a stationary dark soliton (with v% = 0) at a trap center 
and a pair of oscillating solitons (with v\ = —^3) in a cigar-shaped condensate confined in a trap of strength £7 w 0.06. The 
figure shows the evolution of the density as obtained by direct numerical integration of Eq. ©. Additionally, the figure shows 
the trajectory of one of the solitary waves as computed via the equation of motion that makes use of the effective potential of Eq. 
( fTTI i (solid green line) [the other moving dark soliton is the mirror image of the one shown around z = 0, while the third is, by 
symmetry, constrained to stay precisely at z — 0]. It is clear that the agreement between Eq. (0 and the reduced particle picture 
is excellent. 

We complete the analysis of this section by studying the case of dark solitons performing small oscillations around their 
equilibrium positions. In fact, we are going to use the Lagrangian approach devised above to connect the oscillation frequency 
obtained by the multi-soliton dynamics to the eigenfrequencies of the anomalous modes of the stationary soliton states that will 
be obtained by a BdG analysis in the next section. In that respect, it is relevant to consider the simplest case of two well-separated 
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FIG. 3: (Color online) A configuration of three dark solitons, with the central one being stationary and the other two oscillating with a 
frequency u> osc — 0.904Q. Shown are a contour plot depicting the evolution of the density according to Eq. (fJJ illustrating the path of the 
density minimum of one of the oscillating solitons, as well as the respective trajectory calculated by the equation of motion employing the 
effective potential of Eq. ([17} (solid line). The parameter values are N = 2000, uj z = 2n x 53 Hz, uj± — 2nx = 890 Hz (i.e., Q ~ 0.06) 
and initial displacement from the trap center 5zo = 2fim. 



solitons, which are assumed to be almost black (i.e., B\ = B% ~ 1). In such a case, the Lagrangian takes the form, 

r j 2 i „• 2 1 .2 „2 1 .2 „2 n n \ 
L = —Z\ H Z 2 W„„Zi Cl> „Z 9 n ■ (1") 

2 1 2 2 2 osc 1 2 osc 2 sinh 2 ( y /nE(z 2 -z 1 )) 

Then, using the Euler-Lagrange equations, and replacing the hyperbolic function sinh by its exponential asymptote in the case 
under consideration, i.e., for \z 2 — Zi\ >> 0, the following equations of motion are obtained 



Si = -8n 3 /2 e- 2n ° 2 ^-^ -uj 2 oscZi (20) 

y - 8n 3/2 P - 2 "i /2 ( z 2-^l) , ,2 _ 

z 2 — on e - lj osc z 2 . 

The fixed points, Z%, Z 2 of the above system can be easily found by setting the left-hand side equal to zero. The result is: 



Z 2 = -2k = j±=w f^j (21) 



where w(rj) is the Lambert's w function defined as the inverse of 77(10) = we™. Then, considering small deviations (771 , 772) from 
the equilibrium positions Z\, Z 2 , we can Taylor expand the interaction potential keeping only the lowest-order term and, this 
way, derive the following linearized equations of motion: 

m = l6n 2 e-^ z ( V2 - m )-Lv 2 oscm , (22) 

Let us now consider the normal modes of the system and seek solutions of the form T)i = ?/ioe lLJt , i = 1,2, where u> is the 
common oscillation frequency of both dark solitons. Then, substituting this ansatz into Eqs. ( l23b . we rewrite the equations of 
motion as matrix eigenvalue equation, namely: 

W V ~ { 16n 2 e-^ z -lj 2 osc - 16n 2 e~ 4 ^ z J V ' 

To this end, it is possible to obtain from the above system the characteristic frequency u>i = ui osc , which corresponds to in-phase 
oscillations of the two dark solitons, as well as the frequency lu 2 , which corresponds to out-of -phase oscillations of the two 
solitons. The latter is given by: 



lo 2 = ^uj 2 sc + 32n 2 e~^z. (23) 

The above procedure can also be applied to the case of three, almost black, solitons (Bi w 1, i = 1,2,3) considering only 
nearest-neighbour interactions. In this case, the equilibrium positions are given by the following expressions 



1 Von 



Z 2 =0,Z = Z 3 = -Z x = —=w , (24) 
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while the three characteristic frequencies which correspond to the three normal modes of the system are the following, 



In the following section, we will elaborate on the investigation of the stability of the multiple soliton states and the derivation of 
their anomalous modes. 



Having examined the location of the solitary waves and their oscillation eigenmodes and eigenfrequencies in an analytical 
form, we now turn to the numerical investigation of such stationary multi-soliton states, and to their corresponding BdG spec- 
trum. We carry out the relevant computations first for the two-dark-soliton state and subsequently for the three-dark-soliton 
state. 



We start by considering the simplest possible stationary multi-soliton state, namely the second-order nonlinear mode of Eq. 
(fJJ. In the linear limit of N — > 0, this state corresponds to the second-excited state of the quantum harmonic oscillator. The 
excitation spectrum of this state contains a zero eigenvalue (corresponding to the Goldstone mode), two double eigenfrequencies 
located at 57 and 257, as well as infinitely many simple eigenfrequencies. Below we will analyze the excitation spectrum in the 
nonlinear regime (i.e., when the number of atoms N is increased) focusing, in particular, on the two double eigenvalue pairs 
mentioned above. Notice that in our numerical results (see below) we fix u>± to the typical value u±_ — 2tt x 400 Hz. 

First, we consider the two eigenfrequencies located at 57 (for N — 0) which, in the nonlinear regime, obtain opposite Krein 
signature. In particular, one of them has positive Krein signature and corresponds to the dipole mode UauD, while the second 
one has negative Krein signature, i.e., the integral of the norm x energy product, J (\u\ 2 — \v\ 2 )u>dx (in our units), is negative. 
In other words, in the nonlinear regime this eigenvalue becomes the eigenfrequency lq\ of one of the two anomalous modes of 
the system (recall that the number of anomalous modes in the excitation spectrum is the same as the number of dark solitons 
or nodes in the relevant waveform). Notice that both eigenfrequencies originating from 57 in the linear limit [see dashed (solid) 
lines for the one with positive (negative) Krein signature in the top panel of Fig. [4] are real for every value of the number of 
atoms N, indicating the absence of any instability. In fact, the dipolar mode, per the relevant symmetry 10,01, remains fixed at 
u> = 57, while the anomalous mode's frequency in line with the dependence of the single dark soliton mode anomalous mode 
frequency Ulll . 

Next, we consider the eigenvalue pair located at 257 (for N = 0). Similarly to the previous case, these eigenfrequencies 
obtain opposite Krein signature: one eigenvalue has positive Krein signature and corresponds to the background condensate's 
quadrupole mode 10,01, while the second one has negative Krein signature, thus being the eigenfrequency uj2 of the second 
anomalous mode of the system. An important difference from the previous case is that this second pair of double eigenfre- 
quencies does become complex - see bottom panel of Fig. [4] where the imaginary part of these eigenfrequencies is shown 
as a function of N. This implies that the respective nonlinear stationary state is unstable for sufficiently small atom numbers. 
Nevertheless, the instability occurs only near the linear limit (as was also predicted in Ref. fl42ll ). However, as seen in the 
bottom panel of Fig [4] when the number of atoms exceeds a critical value, namely N = 438 for a trap strength 57 = 0.1, all the 
eigenfrequencies become real and the nonlinear state becomes linearly stable. 

We have generally found that the larger the number of atoms, and the stronger the anisotropy of the harmonic trap, the more 
stable the configuration with the two stationary dark solitons is. For example, for 57 = 0.35 this state is unstable up to the number 
of atoms N — 1067, while for 57 w 0.1 the instability occurs for very small condensates, with number of atoms N < 500. Notice 
that a similar behavior was also found in the framework of the ID GP equation considered in Ref. [42] (results not shown here). 
Additionally, regarding the connection of our analysis with experimental observations, we note that within the parameter range 
of relevance to the recent experiment of ll 811 . we found the two-dark soliton state to be linearly stable. 

Let us now return to the excitation spectrum of Fig. |4]and focus on the eigenfrequencies possessing negative Krein signature 
(see solid lines in the top panel of Fig. @), namely the two anomalous modes. In the case of two dark solitons under consideration, 
the physical significance of these two anomalous modes has been discussed in Ref. ll35ll . More specifically, excitation of the 
anomalous mode with the smallest eigenfrequency, u>\, gives rise to an in-phase oscillation, i.e., the two dark solitons move 
towards the same direction without changing their relative spatial separation. On the other hand, excitation of the anomalous 
mode with the largest eigenfrequency, u>2, gives rise to an out-of-phase oscillation, i.e., the two dark solitons move in opposite 




(25) 



IV. STABILITY OF STATIONARY MULTI-SOLITON STATES 



A. The two-dark-soliton state 
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FIG. 4: (Color online) Top panel: the real part uo r of the lowest eigenfrequencies of the second nonlinear mode as a function of the number 
of atoms TV (for a harmonic trap with Q. — 0.1). Eigenfrequencies possessing negative Krein signature, namely eigenfrequencies of the two 
anomalous modes, are depicted by solid lines. Eigenfrequencies possessing positive Krein signature are depicted by dashed lines, with the 
lowest ones starting from u) T = 0.1 and uj r — 0.2 at TV = corresponding, respectively, to the dipole and quadrupole modes. Bottom panel: 
the imaginary part u>i of the eigenfrequency pair starting from ui r = 0.2 at TV = as a function of the number of atoms TV. The vertical dotted 
line at the critical value TV = 438 indicates the splitting of the complex pair into real eigenfrequencies. 



directions with the same velocity and undergo a head-on collision. It should be pointed out that in the context of the analysis 
of the previous section, this information is also encompassed within the eigenvectors r)io of the small perturbations around the 
stationary multi-soliton state. In particular, the two possibilities correspond, respectively, to r)2Q — Via an d V20 — ~Vw- 

The correspondence of the anomalous modes with the normal modes of the two-dark soliton state is confirmed by direct 
numerical simulations. In particular, we have numerically integrated Eq. (O with initial condition the nonlinear stationary 
mode excited by the corresponding anomalous modes, i.e., ip(x; t = 0) = i/jds(x) + u(x) + v*(x). The results (for parameter 
values SI = 0.1 and TV rs 1000) are shown in panels (a), (b) of Fig. In Fig. 0a) (Fig. |2b)), the excitation of the first 
(second) anomalous mode results in an in-phase (out-of-phase) oscillatory motion of the two dark solitons, with the characteristic 
eigenfrequency of the first (second) anomalous mode, i.e., w\ — 0.0784 {u>2 = 0.199). Calculating numerically the maximum 
density of the stationary state and the equilibrium positions of the dark solitons, we find that no — 0.623 and Z = 1.78 ± 0.1. 
This allows us to make a comparison with the ones calculated analytically based on the theoretical approach of the previous 
section. Using Eq. d2.1T >. with u osc = uj\, and Eq. d23l l we find that Z = 1.85 and U2 = 0.205, which differ only 3% from the 
corresponding numerically obtained values. 

Finally, it is worth investigating the manifestation of the instability predicted above for small condensates, due the "collision" 
of the second anomalous mode with the quadrupole mode. In Fig. EJc), we show the evolution of the density of a condensate 
with TV w 400 confined in a trap with strength Q, = 0.1 (as before). For these parameter values, the eigenfrequencies of the two 
anomalous modes are u>i — 0.081 and u>2 = 0.188 4- 0.0025i. The numerical integration of Eq. (O reveals that although the 
initial evolution of the density roughly follows the one observed in Fig. |3b) (up to t w 700), the instability eventually manifests 
itself: the soliton motion excites the quadrupole mode of the system, and this excitation results in a breathing behavior of the 
BEC (see the bottom panel of Fig. [5). 



B. The three-dark-soliton state 

Let us now consider the third-order nonlinear mode of Eq. (|2j. In the linear limit, the excitation spectrum of this state 
consists of a zero eigenvalue, three double eigenfrequencies located at ft, 2£1, and 3f2, as well as infinitely many simple ones. 
In the nonlinear regime, and similarly to the previous case, each of the the three aforementioned pairs obtain opposite Krein 
signature. In Fig. |6]we show the real (top panel) and imaginary (bottom panel) parts of the lowest eigenfrequencies as a function 
of the number of atoms TV (for a trap strength fl = 0.1 as before). In the top panel, the eigenfrequencies depicted by dashed 
(solid) lines correspond to ones with with positive (negative) Krein signature. Regarding the ones with positive Krein signature, 
we note that the lowest ones, starting from uj r = 0.1 and 0.2 for TV = 0, correspond to the dipole and quadrupole modes, 
respectively. As before, the system's ability to sustain dipolar oscillations for all TV with a frequency £7 preserves the dipolar 
frequency at uo r = 0.1 throughout the relevant figure and precludes the possibility of a quartet-inducing collision with the lowest 
in-phase-oscillation anomalous mode of the system. 

We now focus on the two upper double pairs (and their anomalous modes), located at 2S1 and 3f2 in the linear limit. These 
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FIG. 5: (Color online) Spatio-temporal evolution of the condensate density after the excitation of the two anomalous modes. Panels (a),(b) 
show the in-phase and out-of-phase oscillatory motion of the two dark solitons (for S7 = 0.1 and N = 1000); the oscillation frequencies are 
identical to the eigenfrequencies of the first and second anomalous modes, respectively. Panel (c) shows the manifestation of the dynamical 
instability (for SI = 0.1 and N ss 400) due to the collision of the second anomalous mode with the quadrupole mode: the motion of the 
solitons excites the quadrupole mode of the system, resulting in a breathing behavior of the condensate. 
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FIG. 6: (Color online) Top panel: The real part of the lowest eigenfrequencies of the third excited mode (i.e., three-dark-soliton state) as a 
function of the number of atoms N (for a harmonic trap with S7 — 0.1). The notation follows the one used in Fig. [4] The bold circle in the top 
panel indicates the collision of eigenfrequencies of opposite Krein signatures which does not lead to instability (see details in the text). Bottom 
panel: The imaginary parts of the the eigenfrequency pairs starting from ui r — 0.2 (blue) and ui r — 0.3 (green) at N — as a function of the 
number of atoms N. The vertical dotted lines indicate the splitting of the respective complex pairs into real ones. After the second splitting, 
all the eigenfrequencies become real and the three-soliton state is stabilized. 



become complex in the nonlinear regime. As is observed in Fig. |6l the upper pair (starting from 3S1 for N = 0) splits into real 
eigenfrequencies at very small values of the number of atoms N, while the lower pair (starting from 2fi at A~ = 0) remains 
complex for larger values of N. For the assumed trap strength SI = 0.1, the complex eigenfrequencies become real beyond the 
critical value of A" w 880 and, thus, the nonlinear mode becomes dynamically stable. Here, it is worth mentioning that this 
state remains stable for the values of N > 880 considered herein, although at A" w 1395 another collision appears: indeed, 
at a point marked by a circle in the top panel of Fig. [6] the eigenfrequencies starting from u r — 0.3 and 0.4 for N = 0, 
which possess opposite Krein signature, cross each other. Nevertheless, this collision does not lead to instability because the 
eigenmodes associated with these eigenfrequencies remain orthogonal at the collision point. This happens due to the opposite 
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FIG. 7: (Color online) The maximum of the imaginary part of the eigenfrequencies of the third excited state for a fixed number of atoms 
(N — 1000) as a function of the trap strength Q. The insets show the spectral plane for a stable (Q — 0.05) and an unstable (£1 = 0.18) case. 
In the latter case, the complex eigenfrequencies originate from the collision of the second anomalous mode with the quadrupole mode. 

parity of the colliding eigenmodes ll42tl . 

So far, to investigate the stability of the three-dark soliton state (as well as the two-soliton state considered in the previous 
section) we kept the trap strength £1 fixed and varied the number of atoms N. It is also worthwhile (and experimentally relevant) 
to reverse the procedure, i.e., to keep the number of atoms fixed, at N = 1000, and vary the harmonic trap strength £1, as 
shown in Fig. [7] In this figure, it is readily observed that the nonlinear state remains stable up to the critical value £1 = 0.12. 
Beyond this value, the eigenfrequency of the second anomalous mode collides with that of the quadrupole mode and becomes 
complex. Typical examples of the spectral plane, for both stable and unstable cases, are shown in the insets of Fig. [7] For 
the stable case of £1 = 0.05, the eigenfrequencies of the anomalous modes are found to be u>\ = 0.039, 102 = 0.0988, and 
U3 = 0.1626, while for the unstable case of £1 = 0.18 the respective values are ui\ — 0.1489 and L02 — 0.5994; notice that there 
exist also two complex eigenfrequencies at L02 = 0.3379 ± O.Oli, stemming from the above mentioned collision. Once again, 
we calculate numerically the maximum density of the state, no = 0.3817, and the equilibrium positions of the dark solitons: 
Zj = Q,Zs = —Z\ = 4.3±0.1. Using Eq. d24b. with u> n «r = u>i, we find that Z = 4.54, in good agreement with the numerically 
obtained, while the remaining two normal mode frequencies given by Eq. (f26b are found to be U2 = 0.1, 0)3 = 0.1647 and differ 
by less than 2% from the ones obtained by the BdG analysis. 

Similarly to the two-soliton state, we now present results of direct numerical integration of Eq. (O with an initial condition 
given by the third nonlinear mode excited by the corresponding anomalous modes. First we study the stable case, with fl = 0.05 
and N = 1000, and then the unstable case, with £1 = 0.18 and N = 1000. In the stable case, when the first anomalous mode 
is excited, the three dark solitons perform an in-phase oscillation with the characteristic eigenfrequency of the corresponding 
anomalous mode - see Fig. [8l a )- The excitation of the second anomalous mode results in the following configuration: the two 
outer solitons are moving in opposite directions (symmetrically around z = 0), while the center soliton is a stationary one - see 
Fig. Ob). As mentioned in the previous section, this configuration may be experimentally observed following the experimental 
procedure of B18I1 . i.e., by merging two condensates initially prepared in a double-well trap, with a 7r-phase difference between 
them. On the other hand, when the third anomalous mode is excited, the two outer solitons are oscillating in-phase while the 
center soliton is oscillating out-of-phase with respect to the outer ones - see Fig. [SJc)). Finally, in the unstable case, we use as 
an initial condition the third nonlinear state excited by the mode associated with the complex eigenfrequencies. As seen in Fig. 
[3d), initially the outer dark solitons are moving out-of-phase, following the configuration observed in Fig. [8jb). Nevertheless, 
similarly to the two-dark soliton state, the motion of the solitons gradually excites the quadrupole mode of the system, resulting 
in a breathing behaviour of the condensate, signalling the manifestation of the relevant dynamical instability. Clearly, in such a 
case, the oscillation amplitude of the dark solitons is not constant anymore. 

V. EXPERIMENTAL CREATION OF MULTIPLE ATOMIC DARK SOLITONS 

Dark solitons can be created experimentally, e.g., by the method of matter-wave interference, which can be considered as 
a form of density engineering. This method makes use of the fact that an interference pattern in the presence of interatomic 
interactions generates a train of dark solitons B46I1 . Our experimental realization of density engineering involves two BECs, 
initially prepared in a double-well potential. This double well potential is created by the superposition of a crossed optical 
dipole trap and a one-dimensional optical lattice l47ll . Removing the optical lattice leads to the merger of the two initially 
separated BECs in the harmonic trap and to the subsequent creation of an interference pattern which generates the solitons. The 
experimental procedure is described in 11811 . Two oscillating and colliding solitons were observed in that experiment. In the 
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FIG. 8: (Color online) Spatio-temporal evolution of the condensate density after the excitation of the three anomalous modes. In panels (a)-(c) 
stable configurations with Q — 0.05, N = 1000 are shown: in (a) excitation of the first anomalous mode results in a configuration where the 
three dark solitons move together in-phase; in (b) excitation of the second anomalous mode results in a configuration where the outer solitons 
move in opposite directions (each being the mirror image of the other around z — 0) while the middle one is at rest; in (c) excitation of the 
third anomalous mode results in a configuration where the center soliton is moving in opposite direction to the outer dark solitons. Panel (d) 
shows the excitation of the mode associated with a complex eigenfrequency for the unstable case with £1 — 0.18 and N = 1000. Here, a 
dynamical instability manifests itself due to the collision of the second anomalous mode with the quadrupole mode: the motion of the dark 
solitons excites the quadrupole mode of the system, resulting in a breathing behaviour of the entire condensate. 
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FIG. 9: Increasing the distance between the created solitons can be realized by increasing the ramping-down time of the optical lattice tol, as 
shown by numerical simulations of Eq. (O: a) tol = 0ms, b) tol = 2ms, c) tol = 4ms 

following we will extend this scheme to the preparation of a single and of multiple dark solitons in a harmonic trap. 



A. Controlling the created macroscopically excited state 

The fringe spacing of a matter-wave interference pattern depends on the momentum of the two merging atom clouds. There- 
fore, it is possible to vary the number of created solitons in this process by controlling the relative velocity. An increase of the 
relative velocity between the atom clouds lowers the distance between the created solitons and leads to the creation of additional 
solitons further away from the trap center. This can be realized by changing the ramping-down time of the optical lattice tol, as 
shown in Fig. [9] The number and distance of created solitons can also be controlled by the aspect ratio of the trap uj z /uj± or by 
lowering the number of atoms. As thermal fluctuations of the initial phase directly translate into position fluctuations of the dark 
soliton train, care has to be taken to keep this phase as rigid as possible. Since the phase fluctuations scale as oc ksT/Ej 
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(with Ej being the tunneling coupling energy of a double- well system), a careful control of the potential height, i.e. Ej, allows 
to stabilize the phase and thus leads to negligible position fluctuations. This also limits the achievable soliton distances. In our 
experiment we can realize <5$ w 0.06/(2-71"). 

The successful experimental realization of four solitons plus two additional weak ones with extremely high oscillation ampli- 
tude is shown in Fig. [TOt a). The illustrated oscillation dynamics, recorded with a time resolution of 1ms, includes the creation 
process of the solitons starting at the initial double well potential and including the ramping down of the optical lattice. The 
creation process ends at the dashed line, where the final value of the longitudinal trapping frequency is reached after a suitable 
ramping down from the value necessary to obtain the double well potential. 
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FIG. 10: a) Observation of the oscillation of four dark solitons including the creation process. The evolution is averaged over 10 experimental 
runs. The point in time where the creation process of the solitons is finished is marked by the dashed line, b) Experimental observation of three 
dark solitons in a harmonic trap averaged over 16 runs. The creation process of the solitons is not shown in this case. The soliton in the center 
of the trap is at rest (black soliton) whereas the two outer ones oscillate. The time evolution plots were obtained by integrating the images over 
their transverse axis, meaning that each vertical line shows the longitudinal density of the BEC at a certain point in time. 

A matter- wave interference pattern depends on the relative phase difference A<fi of the two merging atom clouds. If A<p = 
a symmetric pattern with an even number of solitons is produced. A small phase difference leads to an asymmetric evolution 
pattern of the created solitons, while a phase difference close to tt leads to the creation of an additional soliton between the other 
ones meaning that an odd number of solitons is produced. Such phase difference can be created by changing the symmetry of the 
potential, which results in an energy difference between the levels of the two wells of the double-well potential. Maintaining this 
asymmetry for a certain hold time accumulates a phase difference between the two BECs. In a simple approximation, the phase 
difference is given by A<f> w AE/h ■ t. By adapting the asymmetry and the time of phase accumulation before releasing the two 
condensates from the double well, arbitrary phase differences can be achieved. Especially interesting is the case were the initial 
phase difference is exactly tt. Then, a black (stationary) soliton is created at the center of the trap, between the oscillating ones. 
Shifting the symmetry of the potential experimentally is realized by shifting the second beam of the dipole trap with respect to 
the optical lattice. 

Fig. [TOl b) shows the experimental realization of three dark solitons in a harmonic trap created by the above discussed method. 
The trap frequencies used in the experiment are v z = 36.1 ± 0.25Hz (longitudinal frequency), v± = 407.5 ± 40.8Hz (transverse 
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frequency). The mean number of atoms in the BEC is N — 1570 ± 146. In this measurement the height of the optical lattice 
is ramped down linearly on a timescale of tql = 2ms. The final value of the longitudinal trapping frequency v z is reached 
after ramping down within 7ms from the value of v™ ltlal — 63Hz necessary for obtaining the double well potential [44]. In the 
performed experiment the oscillation amplitude of the two outer solitons, A osc = (21 ± 0.6)£, is relatively large. Therefore, 
the oscillation frequency is only moderately increased by the soliton-soliton interaction in this case: Vd/v z = 0.775 ± 0.006. A 
numerical simulation yields v c ijv z = 0.761 in good agreement with the experimental result. 

Our method should offer the possibility of creating a single, stationary soliton, corresponding to the first excited state in 
a harmonic trap j45ll . Numerical simulations reveal that this can be achieved by increasing the ramping-down time of the 
optical lattice tol further, which decreases the kinetic energy of the collision process. In the lowest collisional state only one 
interference fringe is produced which produces a single dark soliton. However, for technical reasons [49], the single (moving) 
dark solitons that we created by this method (see Fig. ITTb were fluctuating in position from shot to shot. 

The above results illustrate the possibility of experimentally generating not only a single pair of dark solitons as in fllll . but 
rather of an essentially arbitrary number of such solitons. 
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FIG. 11: (Color online) Single shots of the longitudinal density of the BEC. Increasing the ramping-down time of the intensity of the optical 
lattice from tol = 2 ms to tol = 5 ms offers the possibility of creating a single soliton. Left panel: single shot of the three soliton state 
shown in Fig. [ToJ Note that the stationary soliton in the middle appears with finite density on its minimum in the experimental images. This is 
due to the limited optical resolution. Right panel: single shot of a single dark soliton state. 



VI. CONCLUSIONS 



In the present work, we attempted to quantify the existence, stability and dynamics of multiple atomic dark solitons, by 
examining in detail the prototypical cases of two- and three-dark soliton states. We provided two complemetary viewpoints 
corroborating the same basic picture. A first approach was that of considering the solitons as particles, which interact with each 
other through an exponential tail-tail interaction and are confined within a parabolic trap (appropriately incorporating the effect 
of dimensionality). This particle picture provided us with a detailed understanding of the repulsive nature of the inter-soliton 
interaction and its implication on collision-phenomena and on how it can be combined with the restoring force of the parabolic 
confinement to provide for effective stationary states of the system (i.e., of the dark soliton "crystal"). Within the context of these 
equations of motion the normal modes of this crystal were also examined and were associated with relative motions between 
the solitary waves. The second viewpoint came from the consideration an effective quasi-one-dimensional partial differential 
equation (which incorporates appropriately the transverse confinement of the cloud) and starting from the linear limit of the 
number of atoms N — ► which has the well-known quantum harmonic oscillator eigenfunctions and developing the multi-dark- 
soliton states as natural continuations of appropriate (second- or third- or higher-) excited modes of the linear problem. In that 
context, the excitation spectrum contained the modes of the background BEC (omitted from the particle picture), as well as 
the anomalous modes pertaining to the dark-soliton quasi-particles, which were, in turn, associated with the above mentioned 
normal modes. These two approaches together with a detailed understanding of the experimental setup of lfl8ll provide key 
insights on what types of modes can be excited in the experiment, what intrinsic frequencies should be associated with them 
and, furthermore, what types of instabilities/resonances with background excitations, these modes can be expected to induce. 

One of the future directions of the present program would be to generalize this picture to the extent possible to the n-dark- 
soliton lattice, formulating and addressing questions about the characterization of the normal modes of such a "dark-soliton- 
crystal", as well as questions about the conditions under which this crystal could potentially undergo phase transitions, possibly 
to a state such as a "dark-soliton-gas". On the other hand, another natural generalization of the present program would be that of 
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considering a quasi-two-dimensional analog of the waveforms, the corresponding stability and dynamics, namely that of multi- 
vortex structures and the associated particle picture. Studies along these directions are presently in progress and will be reported 
in future publications. 
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